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1. Introduction 

Quantization of the electromagnetic field in a linear dielectric medium is a nontrivial 
task for various reasons. First of all, since the response of a dielectric to 
external fields is frequency-dependent in general, temporal dispersion should be taken 
into account. The well-known Kramers-Kronig relation implies that dispersion is 
necessarily accompanied by dissipation, so that the quantization procedure has to 
describe an electromagnetic field that is subject to damping. Furthermore, since 
the transverse and the longitudinal parts of the electromagnetic field play a different 
role in the dynamics, the quantization scheme should treat these parts separately. 
For inhomogeneous and spatially dispersive media this leads to complications in the 
quantization procedure, which further increase in the presence of anisotropy. 

When the losses in a specific range of frequencies are small, temporal dispersion 
can be neglected. Field quantization in an inhomogeneous isotropic dielectric medium 
without spatio-temporal dispersion has been accomplished by employing a generalized 
transverse gauge, which depends on the dielectric constant PQ-[B]. 

A phenomenological scheme for field quantization in lossy dielectrics has been 
formulated on the basis of the fluctuation-dissipation theorem [7]-[10]. By adding a 
fluctuating noise term to the Maxwell equations and postulating specific commutation 
relations for the operator associated with the noise, one arrives at a quantization 
procedure that has been quite successful in describing the electromagnetic field in 
lossy dielectrics. An equivalent description in terms of auxiliary fields has been given 
as well [TTJ [12], while a related formalism has been presented recently [13]. However, 
all of these quantization schemes have the drawback that the precise physical nature 
of the noise term is not obvious, since its connection to the basic dynamical variables 
of the system is left unspecified. As a consequence, the status of the commutation 
relations for the noise operator is that of a postulate. 

A justification of the above phenomenological quantization scheme has been 
sought by adopting a suitable model for lossy dielectrics. To that end use has been 
made of an extended version of the Hopfield polariton model [14j in which damping 
effects are accounted for by adding a dynamical coupling to a bath environment. 
Huttner and Barnett [151 116] were the first to employ such a damped-polariton model 
in order to achieve field quantization for a lossy dielectric. Their treatment, which 
is confined to a spatially homogeneous medium, yields an explicit expression for the 
noise term as a linear combination of the canonical variables of the model. In a later 
development, an alternative formulation of the quantization procedure in terms of path 
integrals has been given [T7] , while Laplace transformations have been used to simplify 
the original formalism [TB] . More recently, the effects of spatial inhomogeneities in the 
medium have been incorporated by solving an inhomogeneous version of the damped- 
polariton model [19]-|21j. 

In this way a full understanding of the phenomenological quantization scheme has 
been reached, at least for those dielectrics that can be represented by the damped- 
polariton models mentioned above. The latter proviso implies a limitation in various 
ways. First, one would like to include in a general model not only the effects of 
spatial inhomogeneity, but also those of spatial dispersion. Furthermore, it would 
be desirable to incorporate the consequences of spatial anisotropy, so that the theory 
encompasses crystalline media as well. Finally, while treating temporal dispersion and 
the associated damping, we would like to refrain from introducing a bath environment 
in the Hamiltonian from the start. Instead, we wish to formulate the Hamiltonian in 
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terms of a full set of material variables, from which the dielectric polarization emerges 
by a suitable projection. In this way we will be able to account for any temporal 
dispersion that is compatible with a few fundamental principles like causality and net 
dielectric loss. For a homogeneous isotropic dielectric without spatial dispersion such 
an approach has been suggested before [HI [22] . 

Recently, several attempts have been made to remove some of the limitations 
that are inherent to the earlier treatments. In [23] the effects of spatial dispersion are 
considered in a path-integral formalism for a model that is a generalization of that 
of the original Huttncr-Barnett approach. The discussion is confined to homogeneous 
dielectrics and to leading orders in the wavenumber, so that an analysis of the effects 
of arbitrary spatial dispersion in an inhomogeneous medium is out of reach. In 
|24j crystalline media have been discussed in the framework of a damped-polariton 
model with an anisotropic tensorial bath coupling. A complete diagonalization of the 
model along the lines of [T^l [H] turned out to face difficulties due to the tensorial 
complexity, so that the full dynamics of the model is not presented. Both spatial 
dispersion and anisotropy are incorporated in the quantization scheme discussed in 
|25] . Use is made of a Langevin approach in which a damping term of a specific form 
is introduced. The commutation relations for the noise operator are postulated, as in 
the phenomenological quantization scheme. Finally, several treatments have appeared 
in which a dielectric model is formulated while avoiding the explicit introduction of a 
bath [26] [27] . However, a complete expression for the noise polarization operator in 
terms of the basic dynamical variables of the model is not presented in these papers. 
A direct proof of the algebraic properties of the latter operator is not furnished either. 

In the present paper, we shall show how the damped-polariton model can be 
generalized in such a way that all of the above restrictions are removed. As we 
shall see, our general model describes the quantization and the time evolution of the 
electromagnetic field in an inhomogeneous anisotropic lossy dielectric with arbitrary 
spatio-temporal dispersion. A crucial step in arriving at our goals will be the complete 
diagonalization of the Hamiltonian. It will lead to explicit expressions for the operators 
describing the electromagnetic field and the dielectric polarization, and for the noise 
contribution contained in the latter. In this way the commutation relations for the 
noise operator will be derived rigorously from our general model, instead of being 
postulated along the lines of the phenomenological scheme. Finally, we shall make 
contact with previous treatments by showing how to construct a bath that generates 
damping phenomena in the dynamical evolution of the model. 

2. Model Hamiltonian 

In this section we shall construct the general form of the Hamiltonian for a polariton 
model describing an anisotropic inhomogeneous dispersive dielectric. The result, 
which we shall obtain by starting from a few general principles, will contain several 
coefficients that can be chosen at will. As we shall see in a subsequent section, these 
coefficients can be adjusted in such a way that the susceptibility gets the appropriate 
form for any causal lossy dielectric that we would like to describe. 

The Hamiltonian of the electromagnetic field is taken to have the standard form: 

with the Hermitian vector potential A(r) and its associated Hcrmitian canonical 
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momentum n(r). We use the Coulomb gauge V • A = 0. In this gauge both II 
and A are transverse. The canonical commutation relations read 

[II(r),A(r')] = -i/iST(r-r') , [H(r), II(r')] = , [A(r), A(r')] = (2) 

where the transverse delta function is defined as 6t(t) = I S(r) + VV(47rr) _1 , with I 
the unit tensor. 

The Hamiltonian of the dielectric material medium is supposed to have the general 

form 



H m =h J dr J dujuj C^(r, u) ■ C m (r, 



u) (3) 



with the standard commutation relations for the creation and annihilation operators: 
[C m {T,u),& m (T f ,u/)] = IS{T-T f )5(u-u/) , [C m (r,w),C m (r>')]=0. (4) 

The medium operators commute with the field operators. 

The material creation and annihilation operators are assumed to form a complete 
set describing all material degrees of freedom. Hence, any material dynamical variable, 
for instance the dielectric polarization density, can be expressed in terms of these 
operators. For a linear dielectric medium, the Hcrmitian polarization density is a 
linear combination of the medium operators, which has the general form: 



P(r) = -ifi J dr' dri C m (r',u/) • T(r', r,w') + h.c. 



(5) 



The complex tensorial coefficient T appearing in this expression will be determined 
later on, when the dielectric susceptibility is properly identified. On a par with P we 
define its associated canonical momentum density W, again as a linear combination 
of the medium operators 

W(r) = - J dr' dJ J C m (r', J) ■ S(r', r, J) + h.c. (6) 

with a new complex tensorial coefficient S that is closely related to T, as we shall see 
below. For future convenience we inserted a factor lo' in the integrand and a minus 
sign in front of the integral. 

As W and P are a canonical pair, they must satisfy the standard commutation 
relations 

[W(r),P(r')] = -ifiltf(r-r') , [W(r), W(r')] = , [P(r),P(r')] = 0. (7) 
Hence, the coefficients S and T have to fulfill the requirements: 

dw"f(r",r,w")-T*(r",r / ,o;")-c.c. = (8) 
J dr" J™ duj"cu"S(r",r,uj")-T*(r",r',Lu") + c.c. = \6{r-r') (9) 
dJ'J' 2 S(r",r,w") • S*(r",r', w") - c.c. = (10) 



where the tilde denotes the transpose of a tensor and the asterisk the complex 
conjugate. 

Furthermore, the Hamiltonian should contain terms describing the interaction 
between the field and the medium. Two contributions can be distinguished: a 
transverse part and a longitudinal part. In a minimal-coupling scheme, which we 
shall adopt here, the transverse part is a bilinear expression involving the transverse 
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vector potential A and the canonical momentum density W. To ensure compatibility 
with Maxwell's equations an expression quadratic in A should be present as well, as 
we shall see in the following. For dielectrics with spatial dispersion both expressions 
are non-local. The general form of the transverse contribution to the interaction 
Hamiltonian is 

Hi = -H J dr J dr'W(r)-Fi(r,r')-A(r') + ^ J dr J dr' A(r) • F 2 (r,r') • A(r')(ll) 

with real tensorial coefficients Fi and F2 that will be fixed in due time. In view of the 
form of the second term we may take F2 to be symmetric upon interchanging both its 
spatial variables and its indices. 

The longitudinal contribution to the interaction Hamiltonian is given by the 
electrostatic energy involving the polarization density, which reads 

Here the subscript L denotes the longitudinal part of the polarization density, which is 
defined as [P(r)]^ = —V J cfr'P(r') • V(4 7r |r — r'|) _1 . Furthermore, V is the spatial 
derivative acting on a function of r'. 

The total Hamiltonian H = Hf + H m + Hi + H es is given by the sum of (JTJ) , ([3]), 
(jlip and (fT2"|) . It depends on the tensorial coefficients Fj, F 2 and implicitly on T and 
S through P and W. All of these coefficients can as yet be chosen at will, as long as 
the identities (|8j)— ([TO]) are satisfied. To derive constraints on these coefficients we turn 
to the equations of motion. 

The Heisenberg equations of motion that follow from the total Hamiltonian are 

A(r,t) = -n(r,t) (13) 

£0 

II(r, t) = — AA(r, t) +K [ dr 1 W{r',t) • [Fi(r', r)] T - H [ dr' [F 2 (r, r')] T • A(r', t) 
f-o J J 

(14) 

C m (r,u,t)=-iwC m (T,u,t)-iu J dr' J dr"S*(r,r» • f=x(r',r") • A(r",i) 

+ - f dr'T*(r,r>).[P(r',i)] L , (15) 
£0 J 

where all operators now depend on time. The subscript L' denotes the longitudinal 
part with respect to r'. The time derivative of the polarization density follows by 
combining ([5|) and (fT5j) : 



P(r, t) = -h J dr' J duj' J C m (r', w', t) ■ T(r', r, J) 



(16) 



where h.c. denotes the Hermitian conjugate of all preceding terms. Upon using 
§§§ one finds that the second term (together with its Hermitian conjugate) equals 
— h J dr'Fi(r,r') ■ A(r',t). Furthermore, the last term (again together with its 
Hermitian conjugate) vanishes on account of JHJl. 
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Eliminating II from (|13[) and (114|) we find an inhomogeneous wave equation for 
the vector potential 

AA(r,i) - ^A(r,i) = -fx h J dr'W(r',t) ■ [Fi(r',r)] T 

+ f i nfdr'[F 2 (r ) r')] T -A(r',t) (17) 

where the first term at the right-hand side can be be expressed in terms of the medium 
operators by substituting © . According to the Maxwell equations the vector-potential 
source term, which is given by the right-hand side of (fT7|) , should equal — hq [P(r, t)]x- 
Hence, comparison with (flU]) leads to the identity 

-hjdr'J dr" J duj'u'{C m (r',uj',t)-S(r',r",u}') 
+Ci(r',u',t)-S*(r',r",u')}-[Fi(r",r)] T -h J dr' [F 2 (r, r'] T ■ A(r', t) = 
= -nj dr' J dJ w' {C m (r', u', t) ■ [T(r', r, lo')] t 

+& m (r',J,t) ■[T*(r',r,u>')} T } -h f dr' [F x (r, r')] T • A(r',t). (18) 

Upon equating the coefficient of the vector potential we arrive at the relation 
[Fi(r, r')] TT , = [F2(r, r')} TT , , which connects the transverse parts of the tensors Fi 
and F 2 . A second relation, namely [T(r, r', oj)}t> — J dr" S(r, r", uj) ■ [Fi(r", r')]y , 
follows by equating the coefficient of C m . Combining these two relations with ([8|)- 
(|10[) we thus have found that the tensors T, S, Fi and F 2 have to satisfy five conditions. 
Apart from these constraints the coefficients can be chosen freely while constructing 
our model Hamiltonian. We shall use this freedom to impose instead of the relations 
following from (|18p two somewhat stronger conditions that result upon including the 
longitudinal parts: 

F 1 (r,r') = F 2 (r,r') (19) 

T(r,r»= y < ir"S(r,r", W )-F 1 (r w ,r'). (20) 

As the first of these equalities shows, the bilinear coupling of W with A and the 
quadratic vector-potential term in (|11[) have to occur simultaneously. This is a well- 
known consequence of the minimal-coupling scheme. In view of (|19p we shall omit the 
subscripts of F, in the following. Furthermore, we shall use (|20ll to eliminate S from 
the formalism altogether. 

Summarizing the above results, we have obtained the following Hamiltonian 
for a linear inhomogeneous anisotropic dispersive dielectric interacting with the 
electromagnetic field: 

+ h J dr J dr' J ^ W [C m (r,^)-T(r,r',^) + C ? t „(r, W )-T*(r,r',^)] • A(r') 

+ \h J dr J dr'A(r) • F(r,r') ■ A(r') + ^ J dr {[P(r)] L } 2 . (21) 
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The complex tensorial coefficient T can be chosen freely. It has to satisfy two 
constraints, the first of which has been written already in ©. The second one follows 
by substituting ([2"0j) in (fTUj) : 

j dr" duj" J n f (r", r, J') • T*(r", r', w") - c.c. = 0. (22) 

Finally, insertion of ([2U|1 in @ leads to the equality 

J dv" dio"uj"T{v",v,u")-T*{v" 1 r' 1 uo")+c.c. = F(r,r'). (23) 

This relation defines the real tensor F in terms of T. It shows that F(r, r') satisfies the 
symmetry property F(r, r') = F(r', r), as we know already from the way F2 occurs in 
(jlip . As an integral kernel the tensor F(r,r') is positive-definite. This is established 
by taking the scalar products of (|23p with real vectors v(r) and v(r'), and integrating 
over r and r'. The result is positive for any choice of v. As a consequence, the inverse 
of F is well defined. 

The polarization density is given by ([S]), while the canonical momentum density 
reads according to ([6]) with (|20|) : 



W(r) = -J dr' J dr" J dJ J C m (r', J) • T(r', r", J) ■ V~ x {r" , r) + h.c. (24) 

where the right-hand side contains the inverse of F. 

The Hamiltonian (|2ip has been constructed by starting from general forms for 
its parts Hf, H m , Hi and H cs and requiring consistency with Maxwell's equations. It 
may be related to a Lagrange formalism, as is shown in Appendix A. 

In the following we shall investigate the dynamics of the model defined by (|2~lj) . 
As the Hamiltonian is quadratic in the dynamical variables it is possible to accomplish 
a complete diagonalization. This will be the subject of the next section. 

3. Diagonalization of the Hamiltonian 

We wish to find a diagonal representation of the Hamiltonian ([2~Tj) in the form 

P POO 

H =H J dr J dwuC^r.w) • C(r,u;). (25) 

The creation and annihilation operators satisfy the standard commutation relations 
of the form ([4} . They are linear combinations of the dynamical variables in (f2Tj) : 



(26) 



C(r,w) = j dr'|f 1 (r,r',c;)-A(r')+f2(r,r', W )-n(r') 

P OO 

+ / d W '[f 3 (r,r', W , W ')-C m (r',^)+f4(r,r', W , W ')-CL(r',^)] 
Jo 

with as- yet unknown tensorial coefficients f^, the first two of which are taken to be 
transverse in their second argument. To determine U we use Fano's method |28j : we 
evaluate the commutator [C(r, w), if] and equate the result to Huj C(r, id). Comparing 
the contributions involving the various canonical operators we arrive at the four 
equations 

-f 1 (r,r',c)=cf 2 (r,r',c) (27) 

£0 
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— A'f 2 (r,r',u,)-i/* / dr"f 2 (r,r'».[F(r",r')] T , 
Mo J 

+ J dr" J°° duj"Lo"{f 3 (ry,u;,Lo")-lT*(r",r>,u;")] T , 

-f 4 (r, r", w , w") • [T(r", r', w")] T ,} = w fi(r, r', w ) 
-iW / dr"f 2 (r,r",^)-f(r , ,r", W ')+w'f3(r,r',^,a;') 



|dr"| dr>" J™ du"{f 3 (r,r",u,u;") ■ [T*(r", r'", u")] L „ 



(28) 



h 

+— 

£q J J Jo 

+f 4 (r, r", w , w") • [T(r", r'", W ")M • f (r', r'", u/) = w f 3 (r, r', w , u/) (29) 
-iW j dr"f 2 (r,r",w) • f *(r',r", w') -w'ftfor'.w.w') 

-A | dr "| dv'"l ° O duj"{f 3 (v,r",u J ,w") ■ {T* (v» ,r"> 

+f 4 (r, r", Wj w ") • [T(r", r'", ■ t*(r', r'", «') = w f 4 (r, r', Wj a/). (30) 

The solution of these equations can be obtained by a method that is a 
generalization of that used in our earlier work |19] . The details are given in appendices 
B and C. The results are: 

,2 r 

I l.JI T* I.- .J I \ \f C_// „/ • n\ , 

T> 



h{v,T'^) = ^L J dr"T*(r,r'>).[G(r",r>-iO)b 

f a (r,r / ,w)=i/x w j dr" T*(r,r" ,u>) • [G(r",r', u - i0)] T ' 

f 3 (r,r',cj,a;') = l5(r-r / )<y(w-w')-/«o/iw / dr" / dr'" T(r, r", w) 
■[G(r",r">-iO)] T »' -t(r',r'V) 
+Mo/* ^-TT / dr" [ dr'"T*(r,r'» 

w - u>' - i o y J 

•G(r",r'", W -iO)-f(r',r">') 



(31) 
(32) 



(33) 



f 4 (r,r>,t/)=Mo^ / ^" / dr'" T*(r, r", w) • [G(r", r'", w - i0)] T - • T (r',r'",u/) 



-f^K-^-f / dr" / dr"'T*(r,r", W ) 
w + w J J 



•G(r ,r ,o; - iO) • T (r',r">'). (34) 

The Green function G(r,r',z) occurring in these expressions is defined as the 
solution of the differential equation: 

- G(r, r', z) x VI x V + % G(r, r', z) 
L J c z 

+ J y dr" G(r, r", z) • X (r", r', z) = I 5(r - r') (35) 

The spatial derivative operator V' acts to the left on the argument r' of G(r,r',z). 
According to this inhomogeneous wave equation the Green function determines the 
propagation of waves through a medium that is characterized by a tensor x(r,r',z). 
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The latter plays the role of a non-local anisotropic susceptibility, as will become clear 
in the next section. It is defined in terms of T and its complex conjugate as 

1 



x (r,r\z) = - [ dv" [°° duj 
?o J Jo 

' .f(r",r,w).T(r"y )W ) 



■T(r",r,u,)-T*(r",r>) 



(36) 

LO + Z 

with the frequency argument z either in the upper half or the lower half of the complex 
z-plane, which has got a cut along the real axis. Likewise, the Green function in (|35[) is 
defined in the complex z-plane with a cut along the real axis. Both the susceptibility 
and the Green function are discontinuous across the cut. 

We have succeeded now in finding the diagonal representation of the Hamiltonian 
of our model. The diagonalizing operators are given by P6")) , with coefficients that are 
listed in ([3T |) - (j3"31) . 

4. Field, polarization and susceptibility 

Once we have the diagonal representation of the Hamiltonian at our disposal, we 
can determine the full time evolution of the dynamical variables. In the following 
we will derive the time dependence of the vector potential, the electric field and the 
polarization. As we shall need a few properties of the tensors x an d G, we shall discuss 
these first. 

From its definition (|36|) it follows that the tensor x satisfies the symmetry relations 
x(t,t',z)=x(t',t,-z) (37) 

and 

X*(r,r',z) = X (r,r',-z*) (38) 

so that x(r, r', z) is real on the imaginary axis. The discontinuity across the cut along 
the real axis is given by 

x (r,r> + iO)- X (r,r>-iO) = — /dr"T (r",r,w) • T*(r",r',u;) (39) 

£o J 

for positive u> and by 

X(r,r> + i0)-x(r,r>-i0) = -^^ [ dr"T *(r",r,-w) • T(r",r', -u) (40) 

£o J 



for negative u>. Hence, we may write (|36|) as 

1 f°° 1 

x (v,r',z) = — duj [x(r,r> + iO)-x(r,r', W -iO)j (41) 

2niJ_ 00 uj-z 

which is the well-known Kramers-Kronig relation for the Fourier transform of a 
causal function. The identities ©, (|23[) and can be rewritten in terms of the 
discontinuity across the cut: 

/oo 
dio [ X (r, r', u + iO) - X (r, r', u - iQ)] = (42) 
-OO 

/oo 2 7T i/i 
aio to [ X (r, r', u + i 0) - X (r, r', w - i 0)] = F(r, r') (43) 
-oo £ 

/oo 
dio cu 2 [x(r, r', u + iO) - x(r, r', « - iO)] = 0. (44) 
-oo 
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Incidentally, we remark that for large \z\ the asymptotic behaviour of x follows from 
(gD with (|g5 ]) -([lg j) as 

X (ry,z)^--F(r,r')l + o(lY (45) 
so z \z V 

The above symmetry properties can be used to prove analogous symmetry 
relations for the Green function G. By taking the complex conjugate of ([35]) and 
using ([38]) one derives 

G*(r,r',z) = G(r,r',-z*). (46) 

The adjoint equation of (|35|) reads 

z 2 

- V x [V x G(r, r', z)) + — G(r, r', z) 

<r 

+ J / dv " *( r ' r "' z ) ■ G ( r "> r '' z ) =^{v- r') (47) 

as follows from (|35|) after multiplication by G(r',r",z), integration over r' and a 
partial integration. Comparing this differential equation to that obtained by taking 
the transpose of (|35p and interchanging the position variables one finds with the use 
of (|37f the reciprocity relation: 

G(r,r',*) = G(r',r,-z). (48) 

Having obtained the relevant physical properties of the tensors x an d G, we return 
to a discussion of the time dependence of the dynamical variables. Inverting l|2"6"]) by 
means of the canonical commutation relations we get: 

A(r) =iH J dr' duj T 2 (r', r, w) ■ C(r', w) + h.c. (49) 

Substitution of {32} yields with the help of ((46]) and ([48]) : 



•T(r", r', w) ■ C(r", w) e" 1 " * + h.c. (50) 

where we accounted for the time dependence of the diagonalizing operator C(r",oj). 

By differentiation with respect to space and time we can easily determine the 
time evolution of the electromagnetic fields . Taking the curl of ([50]) we get: 



■f(r", r', u) ■ C(r", uj) &-' lu] 1 + h.c. (51) 



Furthermore, the transverse part of the electric field follows from ([50]) by differentiation 
with respect to t: 



•T(r", r', w) ■ C(r", lu) e~ iuJt + h.c. (52) 

To determine the longitudinal part of the electric field we first have to derive an 
expression for the polarization. 
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The medium operator C TO (r, uj) is a linear combination of the diagonalizing 
operator and its Hermitian conjugate: 

C m (r,w)= J d r'J o [^(r'.r.w'.w) • C(r>') - h(/,T,u/,u) ■ CV.w')] (53) 

as follows by taking the inverse of Substituting and inserting the result 

in (J5J) we get after some algebra 



•i(r"',r",^)-C(r">)e- iwt 

-ihjdr'j dwf(r',r,w)-C(r»e- iwt + h.c. (54) 

where we have employed ((8} and (|36|) . The longitudinal part of this expression can be 
rewritten with the use of (l47l) as: 





Rj£ 









R 







•T(r", r', w) • C(r", w) e~ iuJt + h.c. (55) 

From the Maxwell equation V • (eq E + P) = it follows that the left-hand side 
is proportional to the longitudinal part [E(r, t)]^ of the electric field. The ensuing 
expression for the latter is analogous to (|52l) , so that we arrive at the following result 
for the complete electric field: 

dww 2 G(r,r',w + iO) -f (r",r',w) • C(r",w) e~ iuJt 

+h.c. (56) 

Inspection of (|54|) shows that the polarization consists of two terms. The first 
term is proportional to the electric field, at least in Fourier space and after taking a 
spatial convolution integral. The proportionality factor is x(r, r', u;), which plays the 
role of a susceptibility tensor, as we anticipated in the previous section. The second 
term in |54|) is not related to the electric field. It represents a noise polarization 
density P„(r,<) defined as 

P«(r,t) = -ihjdr'J dwt(r',r,w).C(r',u)e- iwi +h.c. (57) 

that has to be present so as to yield a quantization scheme in which the validity 
of the canonical commutation relations in the presence of dissipation is guaranteed. 
Introducing the Fourier transform P n (r, u>) via 

poo 

P n (r,t)=/ dujP n {r,uj)e~ iuJt + h.c. (58) 
Jo 

and its counterparts E(r, w) and P(r, u>), we get from (|54|) with ([56]): 

P(r,w) = J dr'x(r,r',w + iO) • E(r» + P n (r,w). (59) 

The Fourier-transformed noise polarization density is proportional to the 
diagonalizing operator: 

P n (r, u) = -ih J dr' f (r', r, w) • C(r', w). (60) 
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as follows from ([57)1 and ([55]) . As we have got now an explicit expression for P n (r,w) 
we can derive its commutation relation. By employing (|39p we obtain: 

[P n (r,u,),pt(r>')] [x(r,r> + iO)- X (r,r',a,-iO)] <J( W - W '). (61) 

This commutation relation is a generalization of that postulated in the phenomeno- 
logical quantization scheme for isotropic dielectrics without spatial dispersion [7] -[ID]. 
In the present approach, we have been able to prove its validity. 

Both the fields and the polarization density can be rewritten in terms of P n (r, uj). 
We get from ([ST]) , (|54[) and ([56]) upon eliminating C(r,uj) in favour of P„(r,cj): 

E(r,f) =-n Jdr'J^ duj uj 2 G(r, r', u + iO) • P„(r',w) e~ iw f + h.c. (62) 
B(r,t) =ifj, Jdr'J dwwVxG(r,r> + iO)-P„(r',w)e' iui +li.c. (63) 
P( r)i ) = _i y d r 'y dr"jf dww 2 x(r,r',w + iO)-G(r',r",w + iO)-P n (r",c i ;)e- i " t 
dwP„(r,w)e- iut + h.c. (64) 



By adding ([B"2"]) and (|64[) we get an expression for the dielectric displacement D(r,i). 
Upon using ([47]) we may write it as 



D(r,i) =-Jdr'J dwVx[VxG(r,r> + iO)]-P n (r»e- ik "+h.c. (65) 

Clearly the dielectric displacement is purely transverse. Comparing with (|63|) we find 
that Maxwell's equation V x B(r, t) = /iq <9D(r, t)/dt is satisfied. 

It is instructive to return to the time-dependent representation of the linear 
constitutive relation (E 



P(r, t) = f dr' f dt' X (r, r', t - t') ■ E(r', t') + P„(r, t) (66) 

J J — OO 

with the time-dependent susceptibility tensor defined by writing: 

X(r,r',^ + i0)= / ^x(r,r',t)e iwt . (67) 
Jo 



The convolution integral in the first term of ([66]) . which expresses the causal response of 
the medium, depends on the electric field at all times t' preceding t and at all positions 
r', whereas the second contribution is the noise term, which in classical theory does 
not appear. Sometimes [M] a different splitting of the various contributions to the 
polarization density is proposed, by writing an equation of the general form of (]66[) 
in which the response term covers only a limited range of values of t' , for instance 
t' G [0,t] for t > 0. In such a formulation the convolution integral does not represent 
the full causal response of the medium, so that part of the response is hidden in the 
second term. As a consequence, the latter is no longer a pure noise term, so that it 
cannot be omitted in the classical version of the theory. 

The above expressions for the fields and the polarization density in terms 
of the Fourier-transformed noise polarization density satisfying the commutation 
relations (161|) are the central results in the present formalism for field quantization 
in inhomogeneous anisotropic dielectric media with spatio-temporal dispersion. 
Although we are describing dissipative media, it has not been necessary to explicitly 
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introduce a bath, as is commonly done in the context of damped-polariton treatments 
[T5I IT61 IT91 |2"U] . In the next section, we shall show how a bath may be identified in 
the present model. 



5. Bath degrees of freedom 

In the Hamiltonian (|21| the dielectric medium is described by the operators C m (r,u;) 
and Cj n (r,ct;). The polarization density P(r) and its canonical conjugate W(r) are 
given in ([5]) and (|24p as suitable linear combinations of the medium operators C m and 
C^j. Since the latter depend on the continuous variable u>, they describe many more 
degrees of freedom than P and W. The extra degrees of freedom can be taken together 
to define a so-called 'bath', which is independent of P and W. Although the name 
might suggest otherwise, the bath as introduced in this way is part of the medium 
itself, and not some external environment. Its role is to account for the dissipative 
effects in the dispersive medium, which may arise for instance through a leak of energy 
by heat production. In the following we shall identify the operators associated to the 
bath. Subsequently, we shall show how the Hamiltonian can be rewritten so as to 
give an explicit description of the coupling between the polarization and the bath. In 
this way, we will be able to compare our model to its counterparts in previous papers 

na mi nn hd]- 

The bath will be described by operators Cf,(r, ui) and C|(r, u>) satisfying the usual 
commutation relations. These bath operators are linear combinations of the medium 
operators: 

C 6 (r,w) = J dr' J dJ [H 1 (r,r>, W ')'C m (r', W ') + H2(r,r>, W ')'C|„(r>')] 



with tensor coefficients that will be determined presently. Since the bath variables 
are by definition independent of both P(r') and W(r') for all r', they have to commute 
with the latter. With the use of ([5]) and ([24]) we get from these commutation relations 
the following conditions: 

J dr" l°° du" [Hr(r, r", u, u") ■ T*(r", r V) + H 2 (r, r", u, u") ■ T(r", r', u")] = 

(69) 

J dr" J°° uj" lo" [H x (r, r", u, lo") • T*(r", r',uj") - H 2 (r, r", u, lo") • T(r", r', lo")] = 0. 

(70) 

To determine we start from the following Ansatz: 



H!(r,r>,cy) = J dr" 



8{uj - u') h x (r, r", w ) + \-— h 2 (r, r", u ) 

uj — or + i 



T(r',r",c') 
(71) 

H 2 (r,r', W ,c')-- / dr"— ^h 2 (r,rV).TV,r'V) (72) 
J uj + uj' 

with new tensor coefficients h^. Substituting these expressions in (f6"9"| - (|7T)|) and using 
and (|3^|) . we find that both of these conditions are simultaneously satisfied when 
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hi and h 2 are related as 

J dr" h 2 (r,r'»- X (r",iV + iO) = 

= -L /rfr"h 1 (r,r",c).[ X (r",r',c + iO)- X (r",r',c-iO)]. (73) 

Hence, we are left with a single independent coefficient. It can be determined 
by imposing the standard commutation relation of the form (j4| for C(,(r, u>) and 
c£(r',u/). Using flBSJ) with dTTJ — (JT2J) we arrive at the condition: 

i J dr" J dr'''h 1 (r,r^ W ).[ X (r^r^ W + iO)-x(r'',r''',a;-iO)].h;(r',r''',a;) = 



l(S(r-r'). 



(74) 



The vanishing of the commutator of C&(r, w) with Cf,(r',a/) is warranted on the 
strength of ([73]). 

In view of (l39l) a solution of (1741) is 

(75) 



hi(r,r',w) = f (r'.r.w) 



and hence, on account of ([73]): 

h 2 (r,r» = — [ dr"T*(r,r",Lu) • x _1 (r", r', w + i 0). (76) 

It should be noted that the coefficients are determined up to a unitary 
transformation. This freedom, which is available to as well, corresponds to a 
natural arbitrariness in the choice of the bath operators themselves. 

As the bath operators have been identified now, we can rewrite the Hamiltonian 
so as to clarify their role in the dynamics of our model. To that end we have to 
eliminate the medium operators C m in favor of the bath operators C(,. Employing 
([!]), ([24]) and ([68]) we can write the medium operators as: 



C m (r, w) 
dr' 



dr' T*(r, r', us) 



Y -uj J dr"F- 1 (r',r")-P(r")-W(r') 



diu' 



H;(r',r, W » • C b (r>') - H 2 (r', r, u/, w) • c£(r>') 



(77) 



With the use of this expression the contributions involving the medium operators in 
(|2ip can be rewritten. In this way, we arrive at the following alternative form for the 
Hamiltonian of our model: 

H = / dr {2^ [n(r)]2 + 2^ [VAA(r)]2 } + ^/ dr i dwuC\{r,u>)-C h {r,u) 
+ _££_ J dr J dr' J dr" J dr'"V(r)-?-\r, r') 



F-V.r'") -P(r" 



U dujuj* {x(r', r", to + 10) - x(r', r", w - i 0)]| 
~ J dr {[P(r)] L } 2 + \h J dr J dr' W(r) ■ F(r, r') • W(r') 
h J dr J dv ' w ( r ) ' F ( r > r ') ' A ( r ') + \ h J dv J dr'A{r) ■ F(r, r') • A(r') 
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-'AJdrJ dr' J dr" J dto \c\ (r, lo) • T*(r, r', u) ■ X _1 (r', r", u + i 0) • P(r") 

-h.c] . (78) 

As before, the tensor coefficients T and F are related by ([2"3"]) . while the susceptibility 
X follows from T via ([31>|) . which implies ([3Ti]) . 

We are now in a position to compare the present model with that discussed in 
previous papers [151 HH HH1 HO] • The Hamiltonian ([78]) takes account of anisotropy 
and spatial dispersion. To make contact with the earlier treatments these features 
should be left out. In those circumstances both the susceptibility x an d the tensor 
coefficients T, F are isotropic and local, so that one has for instance: 

X(r,r',z)=x(r,z)\S(r-r'). (79) 

Furthermore, the dielectric medium of the present model has got an arbitrary temporal 
dispersion: the frequency dependence of the susceptibility is governed by that of T, 
as is obvious from (|36jl . In our previous treatment of the inhomogeneous damped- 
polariton model [THl HO] the scalar susceptibility satisfied a sum rule to the effect that 
the integral f Q du> lu 3 [x( r 7 w + i 0) — x(r, u> — iO)], a generalization of which occurs in 
the third term of ([78]) . is proportional to the square of an effective frequency &a(r). 
The latter parameter already figured in the original Hamiltonian in [15] 116] , albeit 
in a space-independent form. Finally, in our earlier work we followed the notation in 
jT3J [T5] by representing the bath operators C(,(r, to) and their Hermitian conjugates 
by equivalent position and momentum operators Y w (r) and Q w (r). Implementing 
this alternative notation here as well, one shows that (|T8[) indeed reduces to the 
Hamiltonian in [19l [20] for the special case of an isotropic spatially-nondispersive 
medium. 



6. Conclusion 



The Hamiltonian model that we have considered in this paper is a suitable tool 
to underpin the quantum formalism for a general linear dielectric medium and 
of the electromagnetic field propagating through such a medium. By solving our 
model we have succeeded in giving a justification of the postulates on which the 
phcnomenological quantization scheme for electrodynamics in dielectric media is 
usually based. 

Our model incorporates many features to warrant the generality of the 
description. Apart from allowing for inhomogeneities and anisotropics of the medium 
it has the virtue of accommodating a quite general spatio-temporal dispersion. In fact, 
the susceptibility tensor of the dielectric medium has been identified in ([56"]) , which 
implies the Kramers-Kronig relation (|4T]) . According to that relation the susceptibility 
in the complex frequency plane is determined by the discontinuity across the cut 
along the real axis. All anisotropic inhomogeneous linear media with spatio-temporal 
dispersion that respond causally to an external electric field are characterized by a 
susceptibility tensor x(r,r',Lj) of the form ([IT]) . Under the assumption that the 
dielectric medium is lossy without net gain, the discontinuity of x(r,r', uj) for u> > 
is a positive-definite integral kernel. Hence, one may use ([39]) to introduce a tensor 
T(r, r',u;), which is uniquely defined up to a unitary transformation. Subsequently, 
one can construct the model Hamiltonian ([2~T]) and proceed with its diagonalization. 
In conclusion, our formalism applies to all anisotropic lossy dielectric media with 
a spatio-temporal dispersion that is compatible with the fundamental principles of 
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causality and positive-definiteness of the dissipative energy loss. Incidentally it may 
be remarked that amplifying dielectric media, which have been treated in the context 
of the phenomenological quantization scheme as well [TUl US], are not covered by 
the present damped-polariton model. To describe media with a sustained gain, e.g. a 
laser above threshold, one has to incorporate a driving mechanism in the Hamiltonian, 
which accounts for the ongoing input of energy that is indispensable for a stationary 
gain. 

As we have shown, the time evolution of the dynamical variables for field and 
matter can be determined completely by deriving the operators that diagonalizc 
the Hamiltonian. The diagonalizing operators are closely related to the noise part 
of the polarization density which plays an important role in the phenomenological 
quantization scheme. The proof of the commutation properties of the noise 
polarization density follows from its relation to the diagonalizing operators. 

In setting up our model Hamiltonian we have avoided to introduce a bath 
environment from the beginning. The subsequent formalism could be developed 
without ever discussing such a bath. Nevertheless, one may be interested in an analysis 
of the complete set of degrees of freedom of the dielectric medium in our model. If that 
analysis is carried out, one finds, as we have seen above, that specific combinations 
of medium variables can be associated to what may be called a bath. The coupling 
of the polarization to this bath can be held responsible for the dissipative losses that 
characterize a dispersive dielectric. 
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Appendix A. Lagrangian formulation 

In this appendix we shall show how the Hamiltonian (f2Tj) can be related to a Lagrange 
formalism. We start by postulating the following Lagrangian for an anisotropic linear 
dielectric with spatio-temporal dispersion that interacts with the electromagnetic field: 



Here A(r) is the transverse vector potential and Q TO (r, u>) are material coordinates 
depending on position and frequency. The polarization density P(r) is taken to be an 
anisotropic and non-local linear combination of these material coordinates of the form 



J Jo 

with a real tensor coefficient To(r, r', u>). One easily verifies that the Lagrangian 
equations have the form 




(A.l) 




(A.2) 



AA(r,t)--A(r,t) 



IH> P(r,t) 



(A.3) 
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Q m (r,uj,t)+ui 2 Q m (r,uj,t) = J dr T (r, r', w) • E(r', t) (A.4) 

with the electric field given as E(r, t) = — A(r, i) — (l/e ) [P(r, i)]_L- The first 
Lagrangian differential equation is consistent with Maxwell's equation, as it should. 
The second Lagrangian equation shows that the material coordinates are harmonic 
variables that are driven by the electric field in an anisotropic and nondocal way. 
Introducing the momenta II(r) and P m (r, u>) associated to A and Q m as 

n(r)=e A(r) (A.5) 

P m (r, u) = Q ro (r, lu) + J dr' T (r, r', w) • A(r') (A.6) 

we obtain the Hamiltonian corresponding to (|A.1[) in the standard fashion. The result 
is: 

^ = / dr {i [n(r)]2 + ^ [VAA(r)]2 } 

/POO 
dr J du {[P m (r, W )] 2 +w 2 [Q m (r,w)] 2 } 

p p p co 

- J dr J dr' J duo P m (r, u) ■ T (r, r', u) ■ A(r') 

+ lJ dr J dv ' J dr " J q ^A(r)-f (r',r,w)-T (r',r",w)-A(r") 

+ ^-|dr{[P(r)] L } 2 . (A.7) 

As a final step we wish to rewrite the Hamiltonian in terms of material creation 
and annihilation operators CJ„ and C m as used in ([2Tj) of the main text. We introduce 
the latter by writing 

P m (r,uO= (jf^j 1/2 Jdr' [C m (r',c).U(r',r,c) + Ct„(r',c).U*(r',r, W )] (A.8) 

1 /2 

Q m (r,w)=i ^) | dr' [C m (r',c).U(r',r,c)-Ct„(r'^).U*(r',r^)] (A.9) 
with tensorial coefficients U that satisfy the unitarity condition 

J dr"U(r",r,cj)-U*(r",r',o;) = I5(r-r'). (A.IO) 
Substituting (|A~8)l and (jAT9|) in (|A~7| we recover (|2ip of the main text, with T given 

by 

T(r,r',w) = -(2Sw)" 1/2 / dr" U(r, r", w) • T (r", r', u). (A.ll) 

Since To is real, one finds that introducing T in this way implies that it fulfills the 
relation 



dr"T(r",r,w) -T*(r",r',w) - c.c. = (A.12) 

which is consistent with (but somewhat stronger than) conditions (JSJ) and (f2"2"|) imposed 
on T in the main text. Adopting the above stronger relation implies that the 
susceptibility (|56")) acquires an additional symmetry property on a par with ([57| and 
(|38p . namely x(r,r',z) = x(r, r', — z). As the validity of (|A.12[) is not essential in 
setting up our Hamilton formalism, we have refrained from using it in the main text. 
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Appendix B. Evaluation of the tensorial coefficients fj 



In this appendix we will show how the equations (|27H30p can be solved. We start by 
using (|27|) to eliminate fi from (|28|) . As a result we obtain the differential equation: 

A'f 2 (r,r',«) + ^-f 3 (r ) r' ) w)-/xo/i f dr"f 2 (r,r",w) • [F(r",r')] T , 

-i/io J dr" du;"u;"{h(r,r",u J ,u;")-[T*(r",r',u J ")} T , 

-f 4 (r,r'>,a/') ■ [T(r",rV)M = 0. (B.l) 

To get an expression for the last two terms we use (|29l) and (l30l) . First, we multiply 
(1291) by T*(r', r"", u>') and integrate over r'. Relabeling the dummy variables we get 



-iW / dr" J dr'" f 2 (r,r" ,uj) -T(r'" ,r," ,lu') -T*(r'" ,r' ,uj') 
-(w-w') J dr"f 3 (ry,to,u;')-T*(r",r',uj') 

+j- J dr" J dr'" J [°° du"{f 3 (T,T",u,w")- [■V*{r",r"',u")] L „, 

+f 4 (r,r", W , W ") • [T(r",r'V')]W ■ T (r"" ( r"', J) ■ T*(r"", r', w ') = 0. (B.2) 

A similar relation is obtained by multiplying by T(r', r"", w') and integrating over 
r': 

- ifi J J dr" J dr'" f 2 (r, r", «) • f * (r'" , r," , w') • T(r"', r', «') 
-(w + w') /" dr"f 4 (r, r",a; )C /)-T(r",rW) 

y W dr'" J dv"" j°° du"{h{vy^^")- [T*(r"yy')] L ,„ 

+f 4 (r,r", W , W ") • [T(r",r"V)]W • f *(r"", r'", c/) • T(r"", r', a/) = 0. (B.3) 
We add (|B.2[) and (|B.3|) . Upon integrating over u' and using © and (f2"3"|) we get 

-i^y dr"f 2 (r,r", W )-F(r",r') 

-co y dr" y °° dw" [f a (r, r", «, w") • T*(r", r', w") + f 4 (r, r", w, w") • T(r", r', J')\ 

+ J dr" y°° dco" lo" [f 3 (r,r^ W ,c,'').T*(r'',r',c,'')-f 4 (r,r'',c,,a;'')-T(r'',r',c,'')] = 0. 

(B.4) 

By taking the transverse part of this relation with respect to r' we get an identity that 
can be used to rewrite (|B.1[) in the form: 



f 2 (r,r',w)x^ x V + ^-f 2 (r,r',w) 

J c z 

-i^ow y dr" f°° duj"{f 3 (r,r",u,u")- [T*(r",r',w")]r' 
+f 4 (r,r", W , W ").[T(r",r', W ")]T'} = 0. 



(B.5) 
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Here we used the transversality of f 2 (r, r', us) in its second argument to write the first 
term as a repeated vector product, with the spatial derivative operator V' acting to 
the left on the argument r' of the function f 2 (r, r', us). 

The integral in (|B.5|) contains the transverse parts of T and T* only. A more 
natural form of the differential equation, with the full tensors T and T*, is obtained 
by introducing instead of f 2 a new tensor g defined as: 

g(r,r',w) ee icjf 2 (r,r',w) 

_1 f dv » [™ du"{h{ry',^u").[T*{v\r>,J>)] L , 
£o J Jo 

+h(r 1 r",Lo,us")-[T(r",r',us")] L ,}. (B.6) 
It satisfies a differential equation, which follows from (|B.5|) as: 

- g(r,r',w)xV x V + -5- g(r, r', us) 
L J c 2 - 

+Vo us 2 J dv" p dus" [f 3 (r, r", us, us") • T*(r", r', us") 

+h(r,r",us,us")-T(r",r',us")] = 0. (B.7) 

The integral contribution still depends on f 3 and f 4 , so that the differential equation 
is not yet in closed form. However, we may rewrite the integral in such a way that its 
relation to g becomes obvious. This can be achieved with the help of the identity: 

J dv" J™ duo" [f 3 (r, r", us, lo") ■ T*(r", r' , us") + f 4 (r, r" , us, oj") ■ T(r", r', J')] = 

= e J dr"g(r,r",«)-x(r",r',w-i0)+8(r,r , ,w). (B.8) 

which contains a tensor s(r,r', ui) that arises while avoiding a pole in the complex 
frequency plane, as we shall see below. Furthermore the right-hand side contains the 
susceptibility tensor x that has been defined in . In (|B.8|) the frequency is chosen 
to be in the lower half of the complex plane just below the real axis. Correspondingly, 
the term — i is an infinitesimally small number on the negative imaginary axis. 

To prove CEl we divide (|B~2| by oj' — lo + iO, with iO an infinitesimally small 
imaginary number. The result is: 

-ifi , U '. n [ dr" f dr"'f 2 (r,r", W ).t(r"',r",c')-T*(r"',r',c') 
uo'-uo + iOJ J 

+ J dr"f 3 (r,r", W , W , )-T*(r"y, W ') 

+— f dr" f dr'" f dr"" f°° duo" {f 3 (r, r", u, lo") ■ [T(r", r '", 

£ lo' - ui + 1 J J J J 



L" 



+f 4 (r, r", lo, co") ■ [T(r", r'", lo")} l „,} ■ f(r"", r'" , lo') ■ T(r"", r', w') 

= 8(co-co')s(r,r',Lo). (B.9) 

At the right-hand side, we have introduced a term proportional to the delta function 
5(u> — u)') to account for the fact that the division by us' — us + i yields a singular 
result for us = us', as discussed in [5H]. The coefficient s is as yet unknown. Likewise, 
upon dividing (|B.3|) by us + us' we obtain: 



us + lo' 



— J dr" J dr"'f 2 (r,r",us).r(r"',r",us')-T(r"',r',Lo') 
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-J dr"f 4 (r, r", w,a/) • T(r",r',a/) 

I rfr" / dr"' / dr"" f°° dco" {f 3 (r, r'>, w") • [T*(r", r">")] £ /« 

+f 4 (r,r",c,c") • [T(r",r"'V)] L ,„} • f r'", u>') ■ T(r"W) = 0. (B.10) 

We subtract (|B.10|) from (|B.9|I and integrate over a/. Inspecting the contributions from 
the last terms at the left-hand sides we find it useful to introduce the susceptibility 
tensor ([36]). Upon employing moreover ([8]), we get as a result of combining (|B.9jl and 
(iBJOl) : 

-ie a uj J dr"f 2 (r,r",^)- X (r",r>-iO) 

+ J dr" J™ eko" [f 3 (r, r", w, u") ■ T*(r", r', J 1 ) + f 4 (r, r'>, w ") ■ T(r", r', 
+ J dv"J dr'"| 00 ^"{f3(r,r", W , W ")-[T*(r",r'", W ")] i „, 

+f 4 (r, r", lu, lo") ■ [T(r", r'", w ") W ' x(r"', r', w - iO) = s(r, r', w ). (B.ll) 

When (|B.6|1 is invoked to eliminate f 2 in favour of g, we recover (|B.8j) . 

Having established (|B.8|) we insert it in (|B.7j) so as to arrive at an inhomogeneous 



wave equation for g: 

- g(r, r', a;) X V x V + — g(r,r',u;) 



+ ^ 1 dr" g(r, r", w ) • x(r", r', w - iO) = - M o ^ 2 s(r, r', «). (B.12) 

To solve g from this differential equation we employ the Green function G(r, r',z) 
associated to the operator at the left-hand side. It has been defined in (|35|) . In terms 
of this Green function the solution of (|B.12jl reads: 

g(r,r',w) = -(J. v 2 J dr"s(r,r",o;)-G(r",r / ,w-iO). (B.13) 

With the use of the above expression for g in terms of s, we can evaluate the 
coefficients fi successively. Let us start with f 2 . From (|B.7|) it follows that the integral 
contribution in (|B.6[) is proportional to the longitudinal part [g(r, r', of g. As 
a consequence, (|B.6[) implies that f 2 equals — {i/ui) [g(r, r', ui)]t'- Hence, upon using 
(|B.13|) and introducing the tensor u by writing s as 

s(r,r',u/) = J dr"u(r,r",£j)-T*(r",r',w) (B.14) 



we obtain: 



f 2 (r,r / ,a;)=i/xow / dr 7 ^ dr'" u(r, r", w) • T*(r", r'", u) ■ [G(r"', r', w - iO)] T '.(B.15) 
On a par with this expression we get from (|27p : 

f^r',^^ y dr" J dr"'u(r,r",c).r(r",r"', W )-[G(r"',r', W -iO)] T '. (B.16) 
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Expressions for f 3 and f 4 follow from (|B.9[) and (|B.10p upon using (|B.13|) (|B.15jl 
and the longitudinal part of (|B.7[) . We get: 

f 3 (r, r', uj, ui') = S(u - ui') u(r, r', ui) 

-^Huo J dr" J dr'" J dr"" u(r, r", ui) • T*(r", r'", u) 
•[G(r"', r"", oj - i 0)] r «« • t(r', r"", u') 

+^ uj ^,_ iQ J dr" J dr'" J dr""u(r,r",u)-T*(r",r'",u) 
•G(r"', r"", w - iO) • t(r', r"", w') (B.17) 
f 4 (r,r / ,w,a/) = Mo S W /" dr"y dr"^ dr""u(r,r",w)-T*(r",r'",cj) 

.[G(r"',r"", W -iO)] T f*(r',r"",c') 

_ Mo £_J^_ J dr" J dr'" J dr""u(r,r",uj)-T*(r",r"',Lo) 

■G^/'.u - iO) • t*(r',r"V). (B.18) 

The introduction of u instead of s has led to the simple form of the first term at the 
right-hand side of (|B.17j) . 



Appendix C. The tensor u 

In appendix B the coefficients fj have been obtained. They are all proportional to the 
tensor u(r,r',o>). In this section, we shall show how this tensor can be determined. 

After insertion of the coefficients in the expression ([26]) it follows that the 
diagonalizing operator C(r, w) itself is proportional to u as well. Since C(r, ui) and 
its Hermitian conjugate must satisfy canonical commutation relations of the general 
form u has to fulfil a constraint that is obtained by evaluating the commutator 
[C(r, ui), CT(r', a/)]. In fact, substituting the expression ([26]) and its Hermitian 
conjugate in the commutator and employing {2j) and ((4|) we arrive at the condition 



dr" 
dr" 



f x (r, r", u>) • ~f 2 (r', r", «') - f 2 (r, r", w ) • f^r', r", a;') 



da/' 



f 3 (r,r", W , W ")-f 3 (r',r",c',c") 
-f 4 (r, r", u, oj") ■ T 4 (r', r", w', = I 5(r — r') 6{u; - J) 



(C.l) 



After the insertion of the formulae for f j we arrive at a bilinear condition for u of the 
form: 

y dr" J dr /// u(r,r'^w)•M(r / ^r // ^a;,a;')•u*(r / ) r /// ,a;') = l^(r-r')^(w-a; / )■ (C.2) 

The explicit form of the tensorial integral kernel M follows by evaluating IjC.ljl 
with (|B.15j) - (|B.18j) . It contains contributions with a variable number of Green 
functions G. The simplest contribution Mo is that without a Green function, which is 
found to be 



M (r,r',w,w') = I <S(r — r') 5{ui — ui 1 ). 



(C.3) 
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The next contribution Mi consists of all terms containing a single Green function: 
M 1 (r,T f ,u,u') = -iMiR J dr" J dr'"T*(r,r'» • jw [G(r",r'V -i0)] r ,„ 



G (r ,r ,a/ - iO) 

, ,'2 



— -G(r",r">-iO) 



w — w' — i 



T" w — u>' — i 
G*(r" / ,r",w / -iO)l-f(r / ,r //, > a/) 



(C.4) 



The first two terms contain the transverse part of the Green function, whereas the 
last two depend on the full Green function. 

Finally, we have got the contributions with two Green functions. Since again both 
the full Green function and its transverse part show up, we can distinguish various 
types of terms. The terms with two transverse Green functions become upon invoking 
©: 

M 2TT {r,r', lu,lu') = ^wu'{u + ui') J dr" J dr'" J dr"" T*(r, r", w) 

• [G(r", r'", w - i 0)] T ,„ • [6* (r"", r>", J - i 0)] ^ • f (r', r"", J). (C.5) 
The terms with one transverse and one full Green function are 
M 2T (r,r',u;,a/) = ^ojJ J dr" J dr 1 " J dr"" J dr v T*(r, r", «) 

• G(r", r'", w - iO) • x(r"', r"", w - i 0) • [g* (r",r"", w' - i 0)] w< 

W [G(r",r"> - i0)] T ,„ • *VW - iO) • G V, r"", W ' - iO)} • f (r',rW) 

(C.6) 



where we introduced the susceptibility tensor (|36|) . The last set of terms we have to 
consider are those with two full Green functions. Again using (|36|) we get 

M a (r,r',a;, tt /) = ^ ft> ff_ iQ j d *" j d *"' J d *"" J ^T(r,r'» 
•G(r",r"> - iO) • [x(r"',r'"> - iO) - xV'", r'", J - iO)] 



•G (r t, ,r"",w / - iO) • T(r , r",a/). 



(C.7) 



Having obtained all contributions to the integral kernel M we are now in a position 
to evaluate their sum. We start by investigating the terms with transverse Green 
functions, as given by (|C.5[) . (|G6[) and part of (|C.4[) . Taking all terms together we 
may write them as 



fiohj dr" J dr"' J dr""T*(r,r'» • j 



-u [G(r",r">-iO)] T „, 



1 5{r"' - r"") - =V G (r"",r"', u' - iO) 



dr v x*(r v , r'" , J - i 0) • G (r"" , r v , w' - i 0) 
I <5(r" - r'") - ^G(r",r">-i0) 
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~ J dr*G(r",r*,«-iO).x(i 



',w-iO) 



G> 



mi in i 

r ,uj 



LO)] J.ffr'.r"",^). 



(C.8) 



The sets of terms that multiply the transverse Green functions are transverse 
themselves, as follows from (|35p . Since the integral of the scalar product of a 
longitudinal and a transverse function vanishes, we can replace each of these transverse 
Green functions by their full Green function counterparts. Subsequently, upon using 
(|35p and performing partial integrations we may rewrite (jC.8[) in the form 



Mo 



H(u + t/) J dr" J dr'" J dr""T*(r,r",w)-{[G(r",r'",cj-iO) X V 

iO)-T(r',r"",a/). 



•G (r 



"" r">' 



(C.9) 



An alternative form for this expression is found upon splitting it in two terms by 
writing the factor (uo + uo') as the difference of uo 2 /(ui — u>' — i 0) and ui' 2 / (u> — uo' — iO). 
Subsequently, we carry out partial integrations in the first term, while we leave the 
second as it stands. Finally, we use (f35|) to eliminate the double spatial derivatives. 
We end up with a set of terms that precisely cancel (|C.7j) and the remainder of (IC.4|) 
(i.e., the terms without the transverse Green functions). 

Collecting the results we find that we are left with (|G3|) . The result for M is thus 
quite simple: 

M(r,r',a;,u/) = lS(r-r')S(uj-uj'). (CIO) 
As a consequence, condition (|C.2j) becomes: 

J dr"u(r,r",cj) •u*(r',r",cj) = \5{r-r'). (C.ll) 

In other words, the tensorial integral kernel u(r, r',u>) must be unitary. For 
convenience we choose from now on: 

u(r,r',w) = IS(r-r') (C.12) 

so that u is independent of the frequency and diagonal in the spatial variables. A 
different choice for u leads to a unitarily equivalent form of the diagonalizing operator 
C(t,w ), as follows from d26j with (|BTT5|> - (|b7T8|) . Upon inserting |GT2|) in (fBT5|) - 
(IB.18j) we finally arrive at expressions (f3"Tj) -([34 |) of the main text. 

As a final check of the expressions for we may verify that the commutator 
[C(r, u>), C(r',a/)] vanishes for all position and frequency arguments. To that end we 
have to check whether the condition 



in 



dr" 
dr" 



fi(r,r'» • f 2 (r',r",c') -f 2 (r,r" )W ) -~h(r' ,r" ,u>') 



dto" 



f 3 (r,r", W , W ")-f 4 (r',r",^, W ") 
f 4 (r,r", W , W ")-f 3 (r',r", W ', W ")" 



(C.13) 



is satisfied. Along similar lines as above one verifies that this is indeed true. 
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